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FOREWORD 

This  program  was  prepared  by  Mr.  John  C.  Sparks  of  the  Structural 
Integrity  Branch  (FBE) , Structural  Mechanics  Division,  Air  Force  Flight 
Dynamics  Laboratory  (AFFDL) , Wright-Patterson  Air  Force  Base,  Ohio. 

The  work  was  accomplished  under  Project  1367,  "Structural  Integrity 
for  Military  and  Aerospace  Vehicles,"  Task  136702,  "Development  and 
Application  of  Structural  Design  Criteria,"  Work  Unit  13670234,  "De- 
velopment of  Airframe  Structural  Design  Load  Prediction  Techniques 
for  Military  Aircraft."  This  report  covers  work  performed  from  1 
April  1976  to  31  August  1976.  Interim  reports  written  were  AFFDL- 
TM-76-71  in  July  1976  and  AFFDL-TM-77-22  in  February  1977.  Appre- 
ciation Is  extended  to  Capt.  Ken  Sims  for  providing  the  wind  tunnel 
data  and  to  Mr.  Theo  Kotwicka  for  assisting  in  chart  preparation. 

This  technical  report  was  released  by  the  author  in  May  1977. 


AFFDL-TR-77-136 

TABLE  OF  CONTENTS 


SECTION  PAGE 

I INTRODUCTION  1 

II  PROGRAM  DESCRIPTION  3 

III  KEYBOARD  DESCRIPTION  6 

IV  SAMPLE  PROBLEM  14 

V CONCLUDING  REMARKS  23 

APPENDICES 

A PROGRAM  VERIFICATION  25 

B "BASIC"  LISTING  OF  SUBSONIC  WING  40 

REFERENCES  48 


AFFDL-TR-77-136 


LIST  OF  ILLUSTRATIONS 


FIGURE  PAGE 


1 

Panel  Numbering  Scheme 

4 

2 

Overlay  Card  for  Subsonic  Wing 

7 

3 

Keyboard  Logic  for  Subsonic  Wing 

8 

4 

WPLOT  for  the  Sample  Problem 

15 

5 

PPLOT  for  a = 5° 

16 

6 

PPLOT  for  a = 1°  and  5°  Deflections  on  Panels  18 

and  24 

20 

7 

Wind  Tunnel  Model  of  the  Aft  Swept  Wing 

26 

8 

C|_  versus  a for  the  Aft  Swept  Wing  at  Mach  = 0.7 

27 

9 

Cq  versus  a for  the  Aft  Swept  Wing  at  Mach  =0.7 

28 

10 

C,  versus  a for  the  Aft  Swept  Wing  at  Mach  =0.8 

29 

11 

Cp  versus  a for  the  Aft  Swept  Wing  at  Mach  = 0.8 

30 

12 

C|^  versus  a for  the  Aft  Swept  Wing  at  Mach  = 0.9 

31 

13 

Cjj  versus  a for  the  Aft  Swept  Wing  at  Mach  = 0.9 

32 

14 

Wind  Tunnel  Model  of  the  Forward  Swept  Wing 

33 

15 

C|^  versus  a for  the  Forward  Swept  Wing  at  Mach  = 

0.7 

34 

16 

Cjj  versus  a for  the  Forward  Sweot  Wing  at  Mach  = 

0.7 

35 

17 

C|^  versus  a for  the  Forward  Swept  Wing  at  Mach  = 

0.8 

36 

18 

Cp  versus  a for  the  Forward  Swept  Wing  at  Mach  = 

0.8 

37 

19 

versus  a for  the  Forward  Swept  Wing  at  Mach  = 

0.9 

38 

20 

Cp  versus  a for  the  Forward  Swept  Wing  at  Mach  * 

0.9 

39 

Vi 


AFFDL-TR-77-136 


SECTION  I 
INTRODUCTION 


A calculator  program  that  analyzes  airloads  on  a wing  is  desirable 
for  quick-look  preliminary  design  work.  The  program  presented  solves 
the  perturbation  velocity  flow  equation 

(i_m2)  = 0 

•dy  32 

for  M < 1 where  (Ji  is  the  velocity  potential.  Wings  are  assumed  to  have 
zero  thickness  and  no  dihedral.  The  program  is  written  in  BASIC  compu- 
ter language  (Reference  1)  and  is  designed  specifically  for  the  Hewlett- 
Packard  9830  or  9831  desk-top  calculator  systems.  SUBSONIC  WING  is  the 
program  name. 


The  starting  point  for  designing  a wing  airloads  prediction  program 
was  to  modify  an  existing  BASIC  routine  that  used  Weissinger's  lifting 
line  theory  (Reference  2).  The  resulting  program  would  perform  an  analy- 
sis in  three  to  four  minutes  but  would  not  consider  changes  due  to  flap 
deflections  or  camber  slope  distributions.  Reference  3 contains  more 
information  on  the  original  Weissinger  program. 


Paneling  techniques  used  in  the  calculation  of  detailed  pressure  dis- 
tributions for  complex  planforms  were  examined  for  possible  adaptation 
to  desk-top  calculators.  An  existing  FORTRAN  routine  (Reference  4) 
developed  by  the  Vought  Corporation,  Systems  Division  (formerly  Ling- 
Tempco-Vought  Corporation  (LTV)),  was  the  baseline  for  the  effort.  This 
program  used  Woodward's  vortex  paneling  method  to  analyze  multiple 
lifting  surfaces  for  both  the  supersonic  and  subsonic  regimes. 
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Two  major  modifications  were  made  to  the  LTV  program  to  decrease 
the  required  computation  time.  The  capability  to  calculate  supersonic 
airloads  was  deleted.  Dihedral  was  deleted.  This  change  eliminated 
certain  coordinate  transformation  logic,  thus  resulting  in  the  assump- 
tion of  only  one  velocity  x:omponent. 

The  final  version  of  the  revised  LTV  program  was  capable  of  analyz- 
ing one  four-point  wing.  This  version  was  translated  into  BASIC  and 
implemented  on  the  HP9830  using  the  Special  Function  Keys  (Reference  5). 
The  Special  Function  Keys  are  a feature  of  the  HP9830  that  eliminates 
the  need  for  programming  flags.  This  is  done  by  allowing  the  programmer 
to  substructure  a code  into  manually  called  subroutines.  SUBSONIC  WING 
was  substructured  into  fourteen  manually  called  subroutines.  The  sub- 
routines were  written  to  maximize  code  operating  efficiency  and  to  mini- 
mize storage  requirements.  Matrix  algebra  was  performed  by  the  Matrix  ROM 
(Read  ^nly  Memory  - Reference  6).  The  capability  to  produce  plots  were 
added  and  a routine  to  consider  break  point  geometry  (Section  II)  was 
introduced. 

SUBSONIC  WING  is  the  result  of  a programming  effort  that  lasted 
from  1 April  1976  to  31  August  1976.  It  is  presently  designed  for  a 
calculator  system  with  8000  words  of  central  memory.  Expansion  would 
be  difficult  since  very  little  memory  is  left  unused.  If  SUBSONIC  WING 
is  adapted  to  larger  systems,  then  additional  panels  may  be  added  by 
simply  changing  DIMENSION  statemc!-ts. 
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SECTION  II 
PROGRAM  DESCRIPTION 

The  SUBSONIC  WING  calculator  program  solves  the  linearized  form  of 
the  perturbation  velocity  flow  equation 


(1  - M^) 


3x^ 


dy^ 


d^(p 

dz^ 


(1) 


for  wings  of  zero  thickness  in  subsonic  flow.  Figure  1 illustrates  the 
coordinate  system.  The  program  will  compute  and  plot  the  wing  pressure 
distribution  as  a function  of  angle  of  attack  and  camber.  External  loads 
are  calculated  by  integrating  the  pressure  distribution  over  the  planfot^m 
area.  SUBSONIC  WING  requires  8K  of  memory  to  operate  and  is  stored  on 
one  3500  word  file.  Appendix  A is  the  program  listing. 

Figure  1 illustrates  a paneling  method  and  panel  numbering  scheme 
for  a 24  panel  wing.  The  wing  of  Figure  1 is  represented  by  quadri- 
lateral panels  with  streamwise  inboard  and  outboard  edges.  SUBSONIC 
WING  allows  a maximum  of  42  panels.  Section  III  presents  a complete 
description  of  the  geometric  input  data  required  to  panel  a wing  and 
Section  IV  presents  a description  on  how  to  panel  a wing.  The  input 
data  may  have  any  consistent  dimensions  of  length. 

The  solution  method  is  based  on  the  subsonic  flow  model  of  F.A.  Woodward. 
A vortex  distribution  on  each  panel  introduces  a constant  load  due  to 
angle  of  attack  and  camber. 
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Figure  1.  Panel  Numbering  Scheme 


AFFDL-TR-77-136 


Compressibility  effects  are  accounted  through  the  Prandtl -Glauert 
transformation 


(2) 


References  7 and  8 contain  more  detailed  descriptions  of  Woodward's 
method  including  derivations  of  constant  strength  source  and  vortex 
panels. 


An  array  of  aerodynamic  influence  coefficients  [A]  is  computed  and 
inverted  by  special  function  key  AERO.  Each  element  a.,  of  [A]  repre- 

■ J 

sents  the  normal  velocity  induced  at  the  control  point  of  panel  i by  a 
singularity  of  unit  strength  on  panel  j.  Given  an  angle  of  attack  and 
camber  slope  distribution,  a system  of  linearized  equations  are  solved 
using  the  flow  tangency  boundary  conditions  for  the  resultant  panel  load. 
The  matrix  formulation  is 


. [(.)  + (C)] 

is  the  column  vector  of  panel  loads  nondimensional ized  by  the 

dynamic  pressure  q.  fa}  is  a column  vector  with  all  entries  equal  to 
the  angle  of  attack  a.  fC)  is  a column  vector  representing  the  slopes 
of  the  mean  camber  line,  fa)  and  fC)  follow  the  same  sign  convention 
and  both  use  degree  measure.  The  positive  direction  is  leading  edge  up. 

Most  of  the  computation  time  is  spent  deriving  and  inverting  [A]. 

This  matrix  should  be  stored  on  ^ data  tape  for  future  use.  The  f i 1 e 
size  needed  is  2P^  where  P is  the  number  of  panels.  Once  [A]'^  is 
computed,  solutions  for  various  angles  of  attack,  control  surface  deflections, 
or  cambers  require  little  additional  time. 
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SECTION  III 
KEYBOARD  DESCRIPTION 

On  the  HP9830  or  9831  calculator  system,  programs  can  be  substructured 
into  manually  called  subroutines  with  the  aid  of  Special  Function  Keys. 
These  subroutines  are  called  when  needed  by  pressing  the  appropriate  key. 
Programming  can  be  very  efficient  since  flags  that  check  the  call  sequenc- 
ing are  eliminated.  Figure  2 is  the  overlay  card  for  SUBSONIC  WING.  An 
overlay  card  is  a template  which  the  user  places  over  the  Function  Keys 
to  identify  subroutine  names  and  positions  on  the  keyboard.  Figure  3 is 
the  subroutine  calling  sequence  showing  branching  options  and  the  relation- 
ship of  each  subroutine  to  the  program. 

A description  of  inputs  and  restrictions  follows  for  each  key. 

SUBSONIC  WING  will  help  guide  the  user  up  to  and  through  the  first  call 
of  routine  DELTACP  by  displaying  the  name  of  the  next  key  to  be  pressed. 


START 

Push  RUN,  then  START  to  initialize  arrays.  "SELECT  CODE"  will  be 
displayed,  indicating  that  the  user  is  to  enter  the  number  of  the  primary 
output  device.  The  second  display  will  then  read  PRESS  SPAN. 

SPAN 

Press  SPAN  to  enter  span  parameters.  Displays  and  descriptions 
are  given  below  in  sequence. 


Display 


Description 


NUMBER  OF  SPAN  STATIONS  The  number  of  stations  needed  to 

define  the  spanwise  panel  bound- 
aries. The  wing  root  is  station 
1 and  the  last  station  is  the 
wing  tip 
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SUBSONIC  WING 


START 

AERO 

1 

SPAN 

CHORD 

break 

MESH 

WPLOT 

AXIS 

PPLOT 

I 

ZERO 

CAMBER 

ALPHA 

DELTACP 

LOADS 

Figure  2.  Overlay  Card  for  SUBSONIC  WING 
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Figure  3.  Keyboard  Logic  for  SUBSONIC  WING 
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Display 


Description 


ENTRY  "I"  Y/S  Span  stations  as  a fraction 

of  the  semi -span.  0 and  1 
will  always  be  the  first  and 
last  inputs.  Inputs  must  be 
in  ascending  order. 


SEMI-SPAN 

PRESS  CHORD 


The  length  of  the  wing  semi- 
span. 

CHORD  is  the  next  key  to  be 
pressed  in  the  input  sequence 
(Figure  3). 


CHORD 


Press  CHORD  to  enter  chord  stations.  Displays  and  descriptions 
are  given  below  tn  sequence. 


Display  Description 

NUMBER  OF  CHORD  STATIONS  The  number  of  stations  needed 

to  define  the  chordwise  panel 
boundaries.  The  leading  edge 
is  station  1 and  the  last  sta- 
tion is  the  trailing  edge. 


ENTRY  "I"  X/C  Chord  stations  as  a fraction 

of  the  local  chord.  0 and  1 
will  always  be  the  first  and 
last  inputs.  Inputs  must  be 
in  ascending  order. 
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BREAK 

Press  BREAK  to  enter  the  break  points.  Break  points  are  the  x 
coordinates  of  leading  and  trailing  edges  for  those  chords  that  define 
a sweep  change.  The  fuselage  centerline  defines  the  positive  x axis 
with  origin  at  the  nose.  Displays  and  descriptions  are  given  below  in 
sequence. 


Display 


Description 


NUMBER  OF  BREAKS  The  number  of  break  point  sets 

(maximum  of  eight)  to  be  entered. 
Two  sets  will  be  needed  to  define 
a four-point  wing. 


SET  "I"  LE,  TE,  YLOC  The  leading  and  trailing  edge 

of  the  chord  at  spanwise  loca- 
tion YLOC.  YLOC  designations 
are  obtained  from  the  table 
generated  by  SPAN.  Break  point 
sets  are  entered  from  inboard 
to  outboard.  The  first  set 
will  always  be  the  coordinates 
of  the  root  chord  leading  and 
trailing  edge  with  a YLOC  value 
of  1.  The  last  set  will  always 
be  the  coordinates  of  the  tip 
chord  leading  and  trailing  edge 
with  YLOC  equal  to  the  number  of 
span  stations. 


PRESS  MESH  MESH  is  the  next  key  to  be  pressed 

in  the  input  sequence  (Figure  3). 


MESH 


Press  MESH  to  generate  the  x and  y coordinates  for  the  panel  corner 
points.  Coordinates  are  generated  from  the  information  entered  in 
SPAN,  CHORD,  and  BREAK.  MESH  must  be  re-pressed  if  the  user  re-presses 
any  of  the  above  keys  to  correct  an  error.  Figure  3 relates  the  logic 
position  of  MESH  to  the  geometry  correction  logic.  The  display  will  read 
PRESS  WPLOT  after  mesh  generation  is  completed. 

10 
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WPLOT 


Press  WPLOT  to  plot  the  paneled  wing.  A HP9862  Plotter  and  a 
HP11271  Plotter  Control  ROM  (Reference  9)  will  be  needed  to  generate 
plots.  The  first  display  will  read  XSCALE,  YSCALE.  XSCALE  should  be 
a value  larger  than  the  absolute  value  of  the  difference  between  any 
two  X coordinates.  YSCALE  snculd  be  a number  larger  than  the  semi-span. 
This  display  will  read  PRESS  AERO  OR  LOADDATA  *,A  after  plotting  is 
completed. 

WPLOT  can  be  used  to  check  and  adjust  wing  geometry.  Figures  1 
in  Section  II  and  Figure  5 in  Section  IV  were  generated  by  WPLOT.  The 
panel  numbers  were  typed  in  manually. 


AERO 

AERO  computes  and  inverts  the  array  of  aerodynamic  influence 
coefficients  [A]. 


The  first  display  reads  MACH  NUMBER.  AERO  accepts  subsonic  Mach 
numbers  only.  Reasonable  results  can  be  expected  for  Mach  numbers  of 
up  to  0.8.  The  second  display  reads  CONTROL  POINT.  This  is  the  down- 
wash  control  point  chosen  as  a fraction  of  the  local  panel  chord  through 
the  centroid.  Reference  7 recomnends  0.95  for  the  control  point.  The 
recommendation  is  based  on  extensive  correlations  of  chordwise  pressure 
distributions  and  lift  slope  curves  for  a variety  of  wing  planforms. 


The  run  time  equation  for  logic  associated  with  the  AERO  sub- 
routine is 

2 3 

T = P + C2  p 


(4) 


where  P is  the  number  of  panels,  T is  the  wall  clock  time  in  seconds, 
and  C.|  and  3'"®  system  constants.  C.|  equals  1.91  and  equals  0.011 

when  running  on  the  HP9830.  On  the  HP9831  , C.j  equals  0.24  and  C^  = 0.001. 
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After  [A]'^  is  generated,  it  should  be  stored  on  tape  for  future 
use.  The  display  will  read  PRESS  STORE  DATA  *,A  THEN  ZERO  and  the 
user  stores  [A]"'  on  file  number  *. 

ZERO 

Press  ZERO  to  set  the  angle  of  attack  and  all  elements  of  the 
camber  array  equal  to  zero.  It  must  be  pressed  after  either  a store- 
data  or  loaddata  command  (Figure  3).  The  key  also  allows  the  user 
to  quickly  correct  cambering  errors.  The  display  will  read  PRESS  ALPHA 
OR  CAMBER  after  the  initialization  is  completed. 

ALPHA 

Press  ALPHA  to  enter  the  wing  angle  of  attack  (ANGLE  OF  ATTACK) 
in  degrees.  The  display  will  read  PRESS  CAMBER  OR  DELTACP  after  the 
entry  is  made. 

CAMBER 

Press  CAMBER  to  enter  panel  deflections  in  degrees.  The  one 
display  reads  PANEL  NUMBER,  CAMBER  which  is  repeated  until  a different 
function  key  is  pressed.  See  Figure  3 for  allowable  options. 

DELTACP 

Press  DELTACP  to  solve  Equation  3 for  the  nondimensional ized 
pressures  ACP/Q.  The  nondimensional ized  pressure  is  printed  for  the 
center  of  each  panel.  The  printing  ord*2r  follows  the  panel  numbering 
scheme  starting  from  the  inboard  leading  edge  and  ending  at  the  out- 
board trailing  edge.  The  display  will  read  PRESS  AXIS  THEN  PPLOT  after 
printing  is  completed. 
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LOADS 

Press  LOADS  to  calculate  wing  loading  from  the  pressure  distribution 

generated  by  DELTACP.  The  external  load,  X center  of  pressure,  average 

chord  length,  and  lift  and  drag  coefficients  are  computed  for  each  bay. 

External  loads  for  the  bays  are  summed  and  a total  is  printed.  Lift  and 

drag  coefficients  are  also  computed  for  the  entire  wing.  The  LOAD/Q 

2 

output  has  the  units  l where  l is  the  unit  of  length.  The  display  will 
read  DONE  after  the  loads  are  computed  and  printed. 

AXIS 

Press  AXIS  to  draw  and  label  the  axes  for  the  chordwise  pressure 
distribution  plot.  The  display  will  read  YMIN,  YMAX,  YTIC  and  the 
user  enters  these  values  to  scale  the  vertical  DELTA-CP/Q  axis.  YTIC 
is  the  marking  increment.  The  second  display  will  read  PRESS  PPLOT 
after  the  axes  are  drawn. 


PPLOT 

Press  PPLOT  to  plot  the  chordwise  pressure  distributions.  The 
one  display  will  read  Y/S,  PRINT  YES{1)  N0(0).  Allowable  Y/S  values 
are  given  in  the  PRESSURE  DISTRIBUTION  printout.  They  will  be  printed 
on  the  plot  for  PRINT  inputs  of  1.  Figures  5 and  6 in  Section  IV  are 
PPLOT  examples.  The  display  is  repeated  until  a different  function  key 
is  pressed  (Figure  3). 
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SECTION  IV 
SAMPLE  PROBLEM 

Program  capability  is  illustrated  in  the  following  problem.  The 
hypothetical  planform  shown  in  Figure  4 is  analyzed  for  the  0.8  Mach 
condition. 

Attention  is  directed  to  several  items  concerning  geometric  input. 
X/C  values  are  progressively  closer  together  as  one  moves  toward  the 
leading  edge.  This  is  done  to  establish  a reasonable  pressure  gradient 
for  force  and  moment  calculations.  The  user  should  avoid  placing  very 
narrow  panels  next  to  wide  panels  when  choosing  Y/S  values.  Panel 
widths  in  the  sample  problem  are  fairly  uniform.  Each  span  station  has 
a YLOC  value  assigned  by  the  program.  Notice  the  YLOC  values  given  in 
the  break  point  entries  correspond  to  Y/S  stations  where  sweep  changes 
occur.  The  planform  is  represented  by  42  panels  which  is  the  maximum 
number  permitted  by  the  program  when  executed  on  an  HP9830  or  HP9831 
with  8K  words  of  central  memory. 

Two  cases  are  examined.  The  first  case  is  for  a = 5“  with  no  con- 
trol surface  deflection.  Figure  5 is  the  PPLOT  for  Case  1.  The  Y/S 
values  printed  on  Figure  5 must  match  those  in  the  DELTACP  printout 
or  an  error  message  results.  The  value  0.55  is  not  printed  which 
illustrates  the  print  suppression  option  in  key  PPLOT.  The  second  case 
is  for  a = 1®  with  5®  deflections  on  panels  18  and  24.  These  two  panels 
define  a trailing  edge  flap.  Figure  6 is  the  PPLOT  for  Case  2.  External 
loads  are  computed  for  both  cases  using  the  LOADS  key.  Notice  the  shift 
in  X center  of  pressure  (X-CP)  caused  by  the  flaps  deflection.  X-CP 
values  are  in  fraction  of  chord  length  measured  from  the  leading  edge. 
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Figure  4.  WPLOT  for  the  Sample  Problem 
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SUBSONIC  WING  PROGRAM 
WING  gp:ometry 


SPAN  STATIONS 
YLOC  Y/S 

1 0 

2 0.15 


3 

4 

5 

6 

7 

8 


0.3 

0.46 

0.64 

0.78 

0.89 

1 


SEI4I-SPAN  * 100 


CHORD  STATIONS 


XLOC 

1 

2 

3 

4 

5 

6 
7 


X/C 

0 

0.05 

0.15 

0.3 

0.5“ 

0.75 

1 


X/C  valves  are  pro- 
gressively closer 
together  as  one  moves 
toward  the  leading  edge. 


BREAK  POINTS 


LE 

TE 

0 

160 

30 

160 

60 

120 

101 

146 

120 

120 

YLOC 

1 

^ (YLOC  entries  correspond 

^ jto  Y/S  values  where 

f 'sweep  changes  occur. 


MACH  NUMBER  = 0.8 
CONTROL  POINT  = 0.95 

RUN  TIME  WILL  BE  4184.208  SECONDS—  HP9830.  The  HP983I 

/run  time  is  515 
I seconds. 
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CAMBER,  ALPHA  SET  EQUAL  TO  ZERO 


ALPHA  - 5 


PRESSURE  DISTRIBUTION 


x/c 

Y/S  = 0.075 

CAMBER 

DELTA-CP/Q 

0.025 

0 

0.492768 

0.1 

0 

0.298366 

0.225 

0 

0.256727 

0.4 

0 

0.229544 

0.625 

0 

0.151585 

0.875 

0 

0.0502885 

X/C 

Y/3  = 0.225 

CAMBER 

DELSa-CP/Q 

0.025 

0 

0.876223 

0.1 

0 

0.478335 

0.225 

0 

0.351449 

0.4 

0 

0.259284 

0.625 

0 

0.146742 

0.875 

0 

0.0580864 

X/C 

Y/S  o 0.38 

CAMBER 

DELTA-CP/Q 

0.025 

0 

1.45718 

0.1 

0 

0.740213 

0.225 

0 

0.46806 

0.4 

0 

0.320125 

0.625 

0 

0.214163 

0.875 

0 

0.119843 

X/C 

Y/S  - 0.55 

CAMBER 

DELTA-CP/Q 

0.025 

0 

1.58344 

0.1 

0 

0.8093 

0.225 

0 

0.50872 

0.4 

0 

0.339638 

0.625 

0 

0.216236 

0.875 

0 

0,108287 
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Y/S  = 0.71 


x/c 

CAMBER 

DELTA -CP /Q 

0.025 

0 

1.67739 

0.1 

0 

0.857417 

0.225 

0 

0.536357 

0.4 

0 

0.342823 

0.625 

0 

0.192704 

0.875 

0 

0.0803978 

Y/S  = 0.835 

X/C 

CAMBER 

DELTA -CP /Q 

0.025 

0 

2.07987 

0.1 

0 

1.04673 

0.225 

0 

0.628225 

0.4 

0 

0.359917 

0.625 

0 

0.166441 

0.875 

0 

0.0632769 

Y/S  « 0.945 

X/C 

CAMBER 

DELTA-CP/Q 

0.025 

0 

3.5046 

0.1 

0 

1.73905 

0.225 

0 

0.99546 

0.4 

0 

0.536499 

0.625 

0 

0.251326 

0.875 

0 

0.104388 

WING  LOADING 
SUMMARY  0?  BAY  LOADS 


Y/S 

I-CP 

CHORD 

LOAD/Q 

CL 

CD 

0.08 

0.26 

145. 

00 

411. ^^6 

0.19 

0.02 

0.23 

0,22 

95. 

00 

352.50 

0.25 

0.02 

0.38 

0.23 

57. 

50 

335.45 

0.36 

0.03 

0.55 

0.22 

52. 

19 

362.10 

0.38 

0.03 

0.71 

0,20 

47. 

19 

255.60 

0.39 

0.03 

0.84 

0.17 

33. 

75 

160.50 

0.43 

O.C.. 

0.95 

0.16 

11. 

25 

85.97 

0.69 

0.06 

TOTAL  LOAD/Q  * 1964,050605 
CL  TOTAL  - 0.295778806 
CD  TOTAL  - 0.025877292 
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CAMBER,  ALPHA  SET  EQUAL  TO  ZERO 


CAMBER  TABLE 
PANEL  CAMBER 

18  5 

24  5 


ALPHA  = 1 


PRESSURE  DISTRIBUTION 


Y/S  - 0.075 


x/c 

CAMBER 

DELTA-CP/Q 

0.025 

0 

0.114711 

0.1 

0 

0.0735524 

0.225 

0 

0.0700804 

0.4 

0 

0.0750052 

0.625 

0 

0.0732862 

0.875 

C 

0.0331229 

Y/S  = 0.225 

X/C 

CAMBER 

DELTA-CP/Q 

0.025 

0 

0.215693 

0.1 

0 

0.122626 

0.225 

0 

0.098886 

0.4 

0 

0.0913831 

0.625 

0 

0.0863851 

0.875 

0 

0.0425873 

Y/S  = 0.38 

X/C 

CAMBER 

DELTA-CP/Q 

0.025 

0 

0.379006 

0.1 

0 

0.202473 

0.225 

0 

0.14376 

0.4 

0 

0.124345 

0.625 

0 

0.138456 

0.875 

5 

0.263446 

Y/S  - 0.55 

X/C 

CAMBER 

DELTA-CP/Q 

0.025 

0 

0.441823 

0.1 

0 

0.239177 

0.225 

0 

0.170909 

0.4 

0 

0.147753 

0.625 

0 

0.160466 

0.875 

5 

0.268631 
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Y/S  -=  0.71 


x/c 

CAMBER 

DELTA-CP/Q 

0.025 

0 

0.48967 

0.1 

0 

0.262799 

0.225 

0 

0.181427 

0.4 

0 

0.138552 

0.625 

0 

0.106183 

0.875 

0 

0.0601908 

Y/S  = 0.835 

X/C 

CAMBER 

DELTA-CP/Q 

0.025 

0 

0.609273 

0.1 

0 

0.314275 

0.225 

0 

0.198272 

0.4 

0 

0.125237 

0.625 

0 

0.0693949 

0.875 

0 

0.0314986 

Y/S  -=  0.945 

X/C 

CAMBER 

DELTA-CP/Q 

0.025 

0 

1.0065 

0.1 

0 

0.502178 

0.225 

0 

0.290816 

0.4 

0 

0.160578 

0.625 

0 

0.0784395 

0.875 

0 

0.0339759 

WING  LOADING 


SUMMARY  or  BAY  LOADS 


Y/S 

X-CP 

CHORD 

LOAD/Q 

CL 

CD 

0.08 

0.33 

145.00 

141.82 

0.07 

0.00 

0.23 

0.31 

95.00 

125.97 

0.09 

0.00 

0.38 

0.42 

57.50 

171.22 

0.19 

0.00 

0.55 

0.40 

52.19 

195.83 

0.21 

0.00 

0.71 

0.26 

47.19 

97.30 

0.15 

0.00 

0.84 

0.20 

33.75 

52.68 

0.14 

0.00 

0.95 

0.17 

11.25 

25.29 

0.20 

0.00 

TOTAL  LOAD/Q  « 810  105971 
CL  TOTAL  - 0.122  ’6547 
CD  TOTAL  - 2.157:  B-05 


shift  caused 
flap  deflection 
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SECTION  V 


CONCLUDING  REMARKS 


The  small  memory  (8-K)  of  the  HP9830  calculator  played  an  important 
role  in  designing  SUBSONIC  WING.  Techniques  were  employed  to  save  core 
not  normally  associated  with  FORTRAN  programming.  Notice  most  of  the 
program  output  is  in  seven  or  eight  digits.  About  90%  of  the  program 
formatting  is  done  using  unformatted  WRITE  statements  coupled  with  TAB 
commands  (Reference  5).  REM  (remark)  statements  are  also  deleted.  Pro- 
gramming done  in  this  fashion  reduced  core  requirements  by  approximately 
15%.  The  core  space  saved  was  used  to  increase  capability  by  enlarging 
storage  and  computational  arrays.  Three  hundred  words  of  core  remain  for 
modifications. 
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APPENDIX  A 
PROGRAM  VERIFICATION 

This  section  compares  program  predictions  to  test  data  obtained 
in  the  Trisonic  Wind  Tunnel  of  the  Aeromechanics  Division  of  the  Air 
Force  Flight  Dynamics  Laboratory.  The  planform  of  the  two  models  tested 
in  the  wind  tunnel  are  depicted  in  Figures  7 and  14.  Figures  8 through 
13  and  Figures  15  through  20  compare  wind  tunnel  results  to  both  USSAEROB 
(Reference  8)  and  SUBSONIC  WING  program  predictions.  The  USSAEROB  model 
included  thickness  effects.  Capt.  Ken  Sims,  of  the  Air  Force  Institute 
of  Technology,  provided  the  tunnel  data  and  the  predictions  from  USSAEROB 
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Figure  8.  C.  versus  a for  the  Aft  Swept  Wing  at  Mach 
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Figure  12.  C.  versus  a for  the  Aft  Swept  Wing  at  Mach 


Figure  15.  C.  versus  a for  the  Forward  Swept  Wing  at  Mach 
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Figure  16.  Cp.  versus  a for  the  Forward  Swept  Wing  at  Mach 
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Figure  17.  C,  versus  a for  the  Forward  Swept  Wing  at  Mach 
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Figure  18.  Cp,  versus  a for  the  Forward  Swept  Wing  at  Mach  = 0.8 
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Figure  20.  versus  a for  the  Forward  Swept  Wing  at  Mach  = 0.9 
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APPENDIX  B 

"BASIC"  LISTING  OF  SUBSONIC  WING 

) 
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START 


10  DIM  ASr42,42],XSr42,4],YSr42,l],ZSr8,3] 

20  Dm  CS[12J,DSL12J,^S[42,1J,PS[42,1J,HSL42,1] 
30  DISP  " SELECT  CODE"; 

40  INPUT  SO 

50  raiTE  (SO, *)TAB25; "SUBSONIC  WING  PROGRAM" 

60  WRITE  (S0,*)TAB80;TAB29;  "V/ING  GEOMETRY "TAB80 
70  DISP  "PRESS  SPAN"; 

80  STOP 


SPAN 


10  DISP  "NUMBER  OF  SPAN  STATIONS"; 

20  INPUT  N1 

30  WRITE  (S0,*)TAB80;TAB29"SPAN  STATI0NS"TAB80;TAB26"YL0C"TAB41 "Y/S" 

40  FOR  1=1  TO  N1 

50  DISP  "ENTRY";!;"  Y/S"; 

60  INPUT  C[I] 

70  WRITE  (S0,*)TAB26;I;TAB41;C[I] 

80  NEXT  I 

90  DISP  "SEMI -SPAN"; 

100  INPUT  Y9 

110  WRITE  (S0,-»^)TAB80; TAB 27 "SEMI -SPAN  =";Y9 
120  DISP  " PRESS  CHORD"; 

130  STOP 


CHORD 


10  DISP  " NUMBER  OF  CHORD  STATIONS  "; 

20  INPUT  N2 

30  WRITE  (S0,*)TAB80; TAB 28 "CHORD  STATI0NS"TAB80;TAB26"XL0C"TAB41 "X/C" 

40  FOR  1=1  TO  N2 

50  DISP  "ENTRY";!;"  X/C"; 

60  INPUT  D[I] 

70  V/RITE  (SO,*)TAB26;I;TAB41  ;D[I] 

80  NEXT  I 

90  DISP  " PRESS  BREAK"; 

100  STCP 
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BREAK 


10  DIS?  " iroi'IBER  OP  BREAKS"; 

20  INPUT  H4 

30  WRITE  (SO,- )TAB80;TAB29 "BREAK  P0INTS"TAB80;TAB22"LE"TAB34"TE"TAB46"yL0C" 

40  FOR  1=1  TO  1T4 

50  DISP  " SET":I;"LE,TE,YLOC"; 

60  INPUT  Z[l,l1,Z[l,2],Z[l-3] 

70  ;7RITE  (S0,*)TAB21 ;ZLI,1J;TAB33;Z[I,2];TAB46;Z[I ,3] 

80  NEXT  I 

90  DISP  " PRESS  MESH"; 

100  STOP 


MESH 


10  N3=(N1-1)*(N2-1) 

20  REDIH  A[N3,N3],E[N3,1],P[N3,1],H[N3,1] 
30  D1=1 

40  FOR  K=1  TO  N4-1 
50  FOR  1=1  TO  N2 

60  Eri,l]=Z[K,l]+D[l>(zrK,2]-Z[K,l1) 

[l,1 J=Z[K+1,l]+DriJ*(Z[K+1,2]-Z[K 


J=Z[K,3]  TO  Z[K+1,3]-1 
=Crz;K+1,3]]-C[Z[K,3]l 
' ^ K,3]] 

TO  N2-1 

=Eri,ii+Fi;i,i 


70  F 

80  NEXT  I 
90  FOR  J= 

100  D2 
110  D3=CLZ 
120  FOR  1= 

130  X[D1 ,1 
140  X D1 ,2] 

150  X'D1,3 
160  X D1,4] 

170  Y Dl ,1 
180  H[D1 ,1 
190  D1=D1+1 
200  NEXT  I 
210  NEXT  J 
220  NEXT  K 

230  DISP  " PRESS  WPLOT  "; 
240  STOP 


K+1,1])-E[I,1] 


=E 

=E 

=E 

=C 

=C 


1,1 j+FLI 

■l+1,l]+Fri+1,ll*( 

;i+i,i]+f!:i+i,i]*( 
J]*Y9 
J+1]*Y9 


J]-D3)/D2 
J+11-D3)/D2 
■crj]-D3)/D2 
C[J+1]-D3)/D2 
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WPLOT 


10  DISP  "XSCAIE  ,YSCALE'' ; 

20  INPUT  S9,S8 

50  SCALE  -0.2*S9.1.2*S9,-0.2*S8,1.2*S8 
40  OFFSET  -Z[1,1J,0 
50  FOR  J=1  TO  N1 
60  D1=INT(J/N1) 

70  PLOT  X[l+(J-1-D1)*(N2-I),1+D1].C[J]*Y9,1 
80  PLOT  X[(J-D1)*(N2-1),3+D1],C[JJ*Y9,2 
90  NEXT  J 
100  FOR  J=1  TO  N2 
110  D1=INT(J/K2) 

120  FOR  K=1  TO  N4-1 

130  PLOT  X[(Z[K,3]-1)*(N2-1)+J-D1 ,1+2*D1] ,01 

140  PLOT  X[(Z[Z+1  ,3]-2)*(N2-1)+J-D1,2*(1+I)1, 

150  NEXT  K 

160  NEXT  J 

170  LABEL  (*,0.9,1.7,0,8/11) 

180  FOR  1=1  TO  N1-2 

190  PLOT  X[1+I*(N2-1),1],C[I+1]*Y9,1 

200  CPLOT  -6,0 

210  LABEL  (*)C[I+1] 

220  NEXT  I 

230  FOR  1=1  TO  N2-2 

240  PLOT  D[I+1]*(Z[1,2]-Z[1,1])+Z[1,1],-0.025*Y9, 

250  CPLOT  -3,0 

260  LABEL  (*)D[l+l] 

270  NEXT  I 

280  DISP  "PRESS  AERO  OR  LOADDATA*  ,A" ; 

290  STOP 


]]*Y9,1 

[K+1,3]]*Y9,2 
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AERO 


10  DISP  "MACH  MJMB3R"; 

20  INPUT  Ml 

50  V/RITE  (S0,*)TAB80;TAB27"MACH  NUMBER  =";M1 
40  M1=SQR(1-Mr2) 

50  DISP  " CONTROL  POINT"; 

60  INPUT  D1 

70  V/RITE  (S0,^)TAB80;TAB26"C0NTR0L  POINT  =";D1 

80  V/RITE  (S0,*)TAB20"RUN  TIME  WILL  BE"  ; (l  .91+0.01 1 *N3)*N3''2 ; "SEC0NDS"TAB80 
90  FOR  1=1  TO  N3 

100  E[I,i1=((1-D1)+(X[I,11+X[I,2])+D1MX[I,3]+X[I,4]))/(2+M1) 

110  F[I,1]=(Y[I,1]+H[I,1])/2 

120  NEXT  I 

130  MT  A=ZER 

140  FOR  1=1  TO  N3 

150  Y1=Y[l,ll 

160  Y2=H[I,1J 

170  FOR  L=1  TO  2 

180  X1=X[I,2^L-1]/M1 

190  X2=X[I,2+L]/M1 

200  D1=(X2-X1)/(Y2-Y1) 

210  S1=SGN(D1)*(1-SGN(D1))+1 
220  D1=AB3(D1) 

230  D2=SQR(1+Dr2) 

240  S2=(-1)'‘L*S1 
250  FOR  J=1  TO  N3 
260  YO=FrJ,ll 
270  X0=e[j,1] 

280  V/0=X0-X1 
290  V/1=X0-X2 

300  A[J,I]=A[J,I]+S2+(FNA(Y0)+PNA(-Y0)) 

310  NEXT  J 
320  NEXT  L 
330  NEXT  I 

340  MAT  A=(Ml/8/PI)*A 
350  MJVT  A=INV(A) 

360  DISP  "PRESS  STOREDATA* ,A  THEN  ZERO"; 

370  STOP 
380  END 

390  DEF  FNA(YO) 

400  Z0=S1*(Y0-Y1} 

410  Z1=S1*(yO-Y2) 

420  D3=SQR(ZO*ZO+WO*WO} 

430  D4=SQR(Z1*Z1+W1*V/1) 

440  D5=(WO+D3)/ZO 
450  D6=(W1+D4)/Z1 

460  K=D2*L0G((D1*W0+ZQ+D2*D3)/(D1*W1+Z1+D2*D4))+D1*L0G(D6*Z0/D5/Z1)+D6-D5 
470  RETURN  K 
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ZERO 


10  KAT  H=ZER 
20  A1=A3=0 

30  raiTE  (S0,*)TAB80;TAB80;TAB20"CAI'mER,  ALPHA  SET  EQUAL  TO  ZERO" 
40  LISP  "PRESS  ALPHA"; 

50  STOP 


CAMBER 

10  V/RITE  (S0,*)TAB80;TAB80;TAB29"CAIIBER  TABLE"TAB80;TAB24"PANEL"TAB40"CAMBER" 
20  DISP  "PAITEL  NUI4BER,  CAMBER"; 

30  INPUT  I,H[I,1] 

40  ’.TRITE  (S0,*)TAB25;I;TAB41;H[I,1] 

50  GOTO  20 
60  STOP 


ALPHA 


10  DISP  "ANGLE  OF  ATTACK"; 

20  INPUT  A1 
30  A3=Al/57. 2957795 

40  V/RITE  (S0,*)TAB80;TAB80;TAB30"ALPHA  ="A1 
50  DISP  "PRESS  CAMBER  OR  DELTACP" ; 

60  STOP 


DELTACP 


10  MAT  E=CON 
20  MiAT  E=(A1)*E 
30  MiAT  P=H+E 

40  MAT  F=( 1/57. 2957795 )*P 
50  MiAT  E=A*F 

60  V/RITE  (SO, *)TAB80;TAB80;TAB24 "PRESSURE  DISTRIBUTION" 

70  D1=1 

80  FOR  J=1  TO  N1-1 

90  VTRITE  (S0,*)TAB80;TAB30"Y/S  ="  ; (C  [J]+C[J+1  ] )/2;TAB80 
100  WRITE  (S0,*)TAB20"X/C"TAB32"CAMBER"TAB42"DELTA-CP/Q" 

110  FOR  1=1  TO  N2-1 

120  WRITE  (S0,*)TAB20;(D[I+1]+D[I])/2;TAB32;H[D1 ,1];TAB42;E[D1 ,1] 
130  D1=D1+1 
140  NEXT  I 
150  NEXT  J 

160  DISP  "PRESS  AXIS  THEN  PPLOT  FOR  PLOT"; 

170  STOP 
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LOADS 


10  WRITE  (S0,*)TAB80;TAB80;TAB29"WING  L0ADING"TAB80 
20  V/RITE  (S0,*)TAB25"SDMMARY  OP  BAY  LOADS" 

30  WRITE  (S0,*)TAB13"Y/S"TAB21"X-CP"TAB28"CH0RD"TAB37"L0AD/Q"TAB48"CL"TAB54''CD" 

40  D1=1 

50  D5=A2=0 

60  POR  K=1  TO  114-1 

70  POR  J=zrK,3]  TO  Z[K+1,3]-1 

80  D2=(C[jJ+C[j+l])/2 

90  DO=(D2-C[z[K,3J])/(C[z[K+1,3]]-C[Z[K,3]]) 

100  B1=M1=0 

110  D3=(Z[K,2]-Z[K,1])*(1-D0)+D0*(Z[K+1,2]-Z[K+1,1]) 

120  POR  1=1  TO  N2-1 

130  B1=B1+E[D1,1]*(D[I+Il-D[l]) 

140  M1=M1+B1*(D[I+1]-D[l]) 

150  D1=D1+1 
160  NEXT  I 

170  D4=Y9*(C[J+1]-C[J])*D3 
180  A2=A2+D4 
190  D6=B1*D4 

200  WRITE  (S0,210)D2;1-M1/B1;D3;D6;B1*C0S(A3);B1*SIN(A3) 

210  PORMAT  13X,P5.2,3X,P5.2,2X,P7.2,P10.2,P6.2,P6.2 
220  D5=D5+D6 
230  NEXT  J 
240  NEXT  K 

250  WRITE  (S0,*)TAB80;TAB22"T0TAL  LOAD/Q  =" ;D5 
260  WRITE  (SO,*)TAB24"CL  TOTAL  =" ;D5/A2*C0S(A3) 

270  WRITE  (SO,*)TAB24"CD  TOTAL  =" ;D5/A2*SIN(A3) ;TAB80 
280  DISP  "DONE"; 

290  STOP 
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AXIS 


10  DISP  "YMIN,YKAX,YTIC"; 

20  INPUT  Y1,Y2,y0 

30  SCALE  -0.3,1.2,1 .2*Y1-0.2*Y2,1.2*Y2-0.2*Y1 
40  YAXIS  -0.1,Y0,yi,Y2 
50  XAXIS  0,1, 0,1 
60  LABEL  (*,1,1.3,0,8/11) 

70  FOR  I=Y1  TO  Y2  STEP  YO 
80  PLOT  -0.1, 1,1 
90  CPLOT  -6,-0. 25 
100  LABEL  (*)I 
110  NEXT  I 

120  PLOT  0.4,-1.5*Y0,1 

130  LABEL  (*,2,1. 5, 0,8/11)  "5a  CHORE" 

140  PLOT  0.15,1.05*Y2,1 

150  LABEL  (*)"CH0RDV/ISE  PRESSURE  DISTRIBUTION" 

160  PLOT  -0.2,0.45*(Y2+Y1)-Y0,1 

170  LABEL  (*, 2, 1.5, PI/2, 8/ll) "DELTA -CP /Q" 

180  LABEL  (*,0.9,1.7,0,8/11) 

190  FOR  1=1  TO  N2 
200  PLOT  Dri],0.1*Y0,1 
210  PLOT  D[IJ,-0.1*Y0,2 
220  CPLOT  -2,-1 
230  LABEL  (*)100*D[l] 

240  NEXT  I 

250  DISP  "PRESS  PPLOT"; 

260  STOP 


PPLOT 


10  DISP  " Y/S, PRINT  YES(1)  N0(0)"; 

20  INPUT  D2,D0 

30  FOR  J=1  TO  N1-1 

40  IF  D2=(C[j]+C[j+l])/2  THEN  60 

50  NEXT  J 

60  D1=(J-1)*(N2-1) 

70  IF  D0=0  THEN  110 

80  PLOT  (D[2]+D[1])/2,E[D1+1,1],-1 

90  CPLOT  -7, -0.25 

100  LABEL  (*)D2 

110  FOR  1=1  TO  N2-1 

120  PLOT  (D[l+l]+D[l])/2,E[D1+I,l],-2 

130  NEXT  I 

140  PEN 

150  GOTO  10 

160  STOP 
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